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Abstract 

Molecular Dynamics simulations of a model bilayer made of surfactant dimers 
in a Lennard- Jones solvent are reported for three sizes of the systems up to an area 
of 100a X lOOcr and for a large interval of the specific areas: from hole formation 
under tension to the floppy state of a buckling compressed bilayer. The transition 
to the floppy state appears quite abrupt and discontinuous; in the floppy state the 
lateral tension is negative. 

Lateral tension and the structure factor were determined for all 3 sizes and all 
areas; the apparent rigidity constant and apparent surface tension are determined 
and correlated with the speciflc area and the flnite size. The replacement of the 
capillary- wave divergence by a pole is accounted for and explained. 
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I. INTRODUCTION. 

The selfassembly of amphiphilic molecules dissolved in a solvent, leads to for- 
mation of bilayers. The properties of these two-dimensional sheets are of great 
interest because of their role in living matter. But it is also of interest to inves- 
tigate the conditions for their formation and equilibrium existence as well as the 
limits of their stability. In particular, the liquid bilayer need not be necessarily 
formed by long-chain molecules; bilayers have been obtained in simulations with 
amphiphilic chains 4 segments long[l,2] and even have been formed by dimers - 
the shortest chain possible [3] . Also, trimers in vacuum have been shown to form 
a stable flat sheet of a bilayer [4]. In [5] we obtained bilayers formed by amphipilic 
dimers, of a new kind, namely "reverse" bilayers, so named by analogy with re- 
verse micelles. Though it may appear that a bilayer can be very easily formed but 
that is not so: in each of these systems the bilayer is formed only within limits of 
temperature, density, and for given interparticle and intermolecular interactions. 
Not much about these limits is known. Always in the same system micelles of 
various shapes can be formed, or the amphiphilic molecules may disperse as a 
solution in the liquid solvent, if the thermodynamic state favors either. Thus the 
limits of existence of bilayers and their stability are of interest. Moreover, we have 
briefly reported[6] for bilayers made of chain molecules unexpected discontinuties 
in the transitions between the extended and floppy states and it is of interest to 
check if such discontinuities also appear in the bilayers made of dimers. Those 
findings [6] for length I = 8 segments were confirmed since for ^ = 4[7]. 

In this paper our previous work on bilayers formed by the shortest chain- 
molecule possible, the dimer[5], is extended to much bigger systems in order to 
study the size dependence. The size dependence study includes the structure 
factor S{q), describing the shape fluctuations of the bilayer. The presence of 
regions of the wave- vector q for which S{q) would have to be negative, and the 
related disappearance of the capillary wave divergence, is explained and resolved. 

In Section II the rather extensive results on the lateral tension are reported 
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and discussed. The transition to floppy buckling bilayer is also examined. The 
example of the apparent discontinuity in the intermolecular energy is given. In 
Section III - the results on the structure factor are given. Section IV is devoted 
to discussion and summary. 

The details of the model and the parameters of the simulations are given in 
Appendix A. 

II. THE LATERAL TENSION AND THE SIZE DEPENDENCE . 

The lateral tension 7 of a bilayer depends on its area A. As the bilayer is 
contained in a square box of volume = L^L^ with periodic boundary conditions, 
the only area available for external control is A = L^Ly = L'^, the edge of the 
box. A is often called the "projected area". The quantity 7 is the free energy 
increase due to area increase at constant volume V, 5F = ^5A. It is computed by 
the same molecular formula as the surface (interfacial) tension between liquid and 
vapor (or any other two fluid phases) - i.e. by the Kirkwood-Buff formula[l-7,8,9] 
as rederived for deformations of a parallelepiped [9]. Constant temperature T and 
particle numbers N, being understood, 



However, the properties of 7 are very much unlike those of the surface tension 
and, partly for that reason, it is called the lateral tension. The specific area per 
surfactant head is defined as 



where is the number of surfactants. Not only does 7 depend on a, but also it 
can take negative values or zero; the state (i.e. the particular area ao = 2Ao/Nd) 
at which 7 = 0, is the "tensionless state". In view of the definition (2.1), at the 
tensionless state 



7 = {dF/dA)v . 



(2.1) 



a = 2A/Nd , 



(2.2) 



dF/dA\A=Ao = 0. 



(2.3) 
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Stability requires (at constant T, V, N again) 



d-f/da > 



(2.4) 



or 



d^F/dA^ > 0. 



(2.5) 



This derivative is related to bilayer lateral (inverse) compressibility [1], Ka, for- 
mally defined as 



It must be positive. But 7 itself can be positive, null, or negative; at positive 
7 the slope of F favors smaller areas (the familiar situation for interfaces with 
surface tension) whereas for negative 7 the bilayer prefers to expand towards the 
tensionless state, (which is the state of lowest free energy, at least locally). 

In our previous simulation work[5] we have produced the values of 7 and thus 
the function 7(a) in series of simulations with different areas, at constant volume 
and at the same thermodynamic state. Unlike in earlier investigations[l,2,3,10], 
we explored [5] negative lateral tensions along with the usual zero and positive 
tensions in a range of areas as large as possible, i.e. often within the entire range 
of stability of the bilayer. The results produced[5] curious and unusual shapes of 
the curves 7(a). Also, in the transition range (near the tensionless area oq) there 
were hints of inflexions in the F{a) curves obtained by numerical integration. It 
was therefore imperative to find which of these new features would be present in 
larger systems, and to study quantitatively the size dependence by extending our 
earlier work to much larger systems. 

Fig. 1 shows the three plots of 7(a) for three systems of nominal sizes Lx ~ 36, 
La; ~ 55, and La; ~ 100 (the unit of length being cr, the collision diameter). Hence 
all areas are in units of cr^. The units used and the parameters are given in all 
detail in Appendix A. 

Remarkably, positive 7's fall on a common curve; there is no visible effect 
of the size. For small areas another region is obtained with a very small, almost 
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Ka = a{d-f/da) > 



(2.6) 



negligible, positive slope and a negative 7. In the latter region a clear size de- 
pendence can be seen: smaller I7I for larger system. The density profiles, i.e. 
the 2;— dependence of concentrations, and the direct imaging of the system spatial 
configuration, show a kind of a rough, fioppy, fuzzy, crumpled, foamy, or buckled, 
state of the bilayer in this region. Often the bilayer is full of solvent particles 
which suddenly find it possible to penetrate into the fioppy bilayer. 

The conclusion is that there are two regions, (E) and (C). In (E), the region of 
large specific areas a, the bilayer is extended and gently undulating, the system is 
stable, 7 > 0, the slope drijda > is positive, large, and common to all sizes. No 
visible size dependence in 7(a) is seen as the data points fall on a common curve. 
In (C), the region of small specific areas a, the bilayer is compressed and fioppy, the 
lateral tension is negative 7 < 0. The derivative d7/da > throughout the entire 
range. This implies stability. In few instances, at the very lowest investigated 
values of a, the sign of the derivative was uncertain in the largest system, still 
practically vanishing within the scatter of the data points. 

At each size there is a smooth transition between (E) and (C) , as can be seen 
in Fig.l. This transitional region is rather narrow. 

The scaling with size of negative tension data in region (C), was found sur- 
prisingly unambiguous. 

The area dependence in (C) is very weak; a plateau value 7p of 7 can be 
assigned to each size. When plotted against the system size, a very satisfactory 
scaling is obtained in Fig. 2. For the measure of the size one can take the edge 
Lx = Ly, or its square, the area A = L^Ly, or Na = 2 A/ a. Clearly 7^ tends to 
zero. Moreover it does so along a straight line 7^ = —const. /A. 

There is one calculation [11] which explicitely predicts a negative lateral ten- 
sion 7 and its saturation to the lower bound given by 

«din + 7 = 0. (2.7) 
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Taking for q^i^ ^^e fowest possible Fourier vector q^^^ — i'^'^/Lx)^ we obtain 

^i f,_ = -K ■ 4:'n'^/A = -87r^Av/iVdO (2.8) 

For an individual system A^^ =const. and 7 ~ const. /a\ otherwise the size de- 
pendence at each specific area o is ~ l/Nd- This is very satisfactory as all these 
predictions of the theory are confirmed by our data. The quantity k, is the theo- 
retical microscopic bending (rigidity) coefficient of the membrane, devoid of any 
surface tension in this theoretical picture. The scaling (2.8) was obtained inde- 
pendently by Otter[7]. 

There are other quantities which show a similar kind of discontinuous jump. 
Most obviously, the lateral inverse compressibility Ka jumps from a high value 
in region (E) - proportional to a as the slope is common to all sizes and constant 
within (E) (see Fig.l) - through a smooth transition to a very low value, still 
positive in the region (C) of the fioppy bilayer. 

The intermolecular energy is another example. Fig.3 shows just one plot of 
U vs. a for the largest system; The plot is in accordance with those of 7(a) in 
Fig.l: two regions of (almost linear) variation of U{a) are joined by a transition. 
The region (C) of negative 7 corresponds to the region of negative slope dU /da 
here and the region (E) of positive 7 corresponds to the region of positive slope 
dU j da. When comparing different sizes, no trend was detected in the slopes nor 
in their difference. The latter appeared to be independent of the bilayer size. 

Also histograms of the angle between the axis of the surfactant dimer and 
the 2;— axis vary with a with some (rounded) discontinuity. Further examples are 
found with the parameters of the structure factor S (q) in section III. 

In small systems the tensionless state 7(00) = belongs to the region (E) 
and the curve 7(a) crosses the ordinate axis with little change in slope. Also Ka 
varies little and its value at the tensionless state can be used for prediction of 7(a) 
at a range of values a > oq. It is not so for large systems, where the tensionless 
state falls within the transition region and the derivative d'y /da changes very fast 
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in the immediate neighbourhood of the tensionless state. This can be seen again 
in Fig.l. 

The apparent transition merits investigation. The derivative of the lateral 
tension w.r.to the area, 7' = dj{a)/da, undergoes a discontinuous jump, just like 
an order parameter would in a first-order transition. This jump is apparent from 
Fig.l and 2 and can also be seen from plots of the computed 7' from raw data 
points. 

However, for a first-order transition the rounding - as resulting from the finite 
size - should disappear with size increasing indefinitely. The area ao of the ten- 
sionless state should stabilize to a definite limit. The point of maximum curvature 
(maximum change in slope) should merge with the tensionless point gq. 

Fig.4 shows the hypothetical diagram of 7 vs. a at constant T, V, N in which 
rounding disappeared. The plot reduces then to two straight lines meeting at a 
transition point (at,7t). A straight line 7(a) = +s{a — ao) with positive slope 
s for a > ao (region E) is continued down to the transition point 7^,0^ with 
ao > atjjt < 0. The other line may be a constant or may follow the prediction of 
eq.(2.8) as —e/a with e— > 0+. This plot would be correct for a classical first-order 
transition, sharp in the limit of macroscopic system. 

The derivative 7 ' may be described by an interpolation formula 

/(o) = (/i + /2)/2 + (/2 - /i) tanh(c (a - a*))/2 . (2.9) 

Here /i, /2,c, a* are parameters and c describes the sharpness of the transition. 
Integration of /(a) produces a function for fitting 7(a) 

9{a) = 90 + (/i + /2)(a - a*)/2 + (^ - /i) log[2cosh(c(a - a*)]/(2c) (2.10) 

The function g{a) produced very good least-square fits of 7(a) for all three sizes 
with c and a* as free parameters. 

However, the sharpness parameter c did not vary significantly with size. The 
other parameters varied very little with size either. Such patterns are not expected 
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for finite-size scafing of first-order transitions for which the sharpness parameter 
increases monotonously with size. 

It therefore appears that (2.9) and (2.10) are just interpolation formulae for 
each individual system. 

The formulae of Reference[ll] which successfully predicted the saturation of 
7 to negative values in the floppy regime of the bilayer, may also be used for 
the entire region of specific areas. The smoothed mesoscopic area (as opposed to 
theoretical intrinsic area A of the membrane) cannot be different from A = LxLy 
as our bilayer is enclosed in a box with periodic boundary conditions. With this 
simplification, but without any change of physical meaning of all quantities, we 
can reproduce not too badly the curves 7(a) as shown in Fig.l. In nondimensional 
form, appropriate for discussion of size effects, eq(14) of Reference[ll] reads: 

U-1 + ua/v2 = (1/c) log[^4^^] 
with the abbreviations 

u = a/a ; c = 87i(3k ; b = Sn'^/dNd . 

It is an equation to be solved for a = a{a) or a = a{a). The parameter l/v2 > 
makes the fixed intrinsic area a slightly elastic. There is no doubt that a of [11] 
is our 7 as it is coupled to A = L^. 

The relevant parameters have reasonable values: k ~ 5 — 10 ~ (2.5 — 5)kT, 
a = 2A/Nd = afixed ~ 1-13, Qmax ~ 27t/w with the width w =4.5, 5, or 6, qmin as 
explained above q^^n = Sir^/Nda, the elasticity parameter is large V2 ~ 30 — 40.. 

However, the agreement is only rough; the transition predicted theoretically 
is much too smooth. The excellent fit of the fioppy region (C) provides reliable 
values of the bending constant n but then a straight line in region (E) is impossible 
to obtain with small values of k. Conversely, large values of k, implying a rigid 
membrane, can produce a straight line in region (E) but then agreement in region 
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(C) is lost. The best overall reproduction of data was obtained for the smallest 
system. Also, we ran into numerical contradictions at the tensionless state; there 
the equations predict a steady shift of oq to lower values with the slope at the 
tensionless state decreasing to zero and our data do not fit this pattern. 

Hopefully, the theory sketched in Reference [11] can perhaps serve as a starting 
point for improvement and also for a prediction of the structure factor. 

It is known[12,13] that many transitions leave a signature on Cy. The heat 
capacity Cy and the function Cy{a) were also extracted from our data. We find a 
constant Cy within the scatter of data, independent of a, for all three sizes. 

III. THE STRUCTURE FACTOR AND ITS POLES . 
A. Extracting and fitting S{q). 

The shape fluctuations of the two dimensional bilayer sheet immersed in three 
dimensions, are similar to capillary waves and it is expected that the capillary 
wave theory may be applicable. However, at the tension free state the lateral 
tension 7 vanishes. It was demonstrated[l], in a simulation experiment, that the 
capillary wave divergence 

Siq) - kT/^q^ (3.1) 
is replaced by a stronger divergence 

S{q) ~ kT/Kq^ (3.2) 

ruled by the rigidity coefficient n. This important result [1] was obtained for the 
tension-free state 7 = 0. In our work[5], we have investigated states with non-zero 
7 and have represented the divergent term as 

S{q) r^l/{kx'^ + gx) x = q^ (3.3) 

The structure factor[l,5,7,13,14,15] for the bilayer is related by Fourier transfor- 
mation to the correlation < h{x,y)h{x' ,y') >. 
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The function h{x, y) describes the shape of the bilayer by giving the "height" 
or the 2;— coordinate, as a function of the position in the a;, y plane. The latter is 
the plane of a perfectly flat bilayer. Now one must define what is meant by the 
position h. To avoid a certain degree of arbitrariness involved in any smoothing, 
we use the actual positions r of each head of the amphiphilic molecule, as 

Therefore in 

5(q) =< /iq/i_q > (3.4) 

we take 

/iq = {l/Na) exp[^qRj] x z^. (3.5) 
i 

Here Rj = [xj.yj) is a two-dimensional position vector of the dimer j, Zj its 
"height", the 2— coordinate, q = (qx^Qy) is the two-dimensional Fourier vector. 
Care is taken of the translational invariance and periodic boundary conditions. 

Unlike here (see also [5,6,7] ), in some simulation work on bilayers the x,y 
grids were introduced together with a recipe by which the local height h{xn, yn) 
would be computed[l,14,15]. 

The average < ... > in (3.4) is the time average over a Molecular Dynamics 

run. 

We found earlier[5] and we find now that the divergent term (3.3) must he 
supplemented by terms describing the smooth and regular background. There are 
bulk contributions which extend down to q— > 0. 

The structure factor S{q) determined from the simulation was fitted to the 
semiempirical formula 

S = 1/ (kx^ + gx) + p/x + wq + wix + W2X^ (3.6) 

with p — or not and with other constants sometimes put to zero. We introduce 
the term p/x after Reference[l] where the form 1/kx^ +p/x was used. Our results 
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for S{q) were obtained afong with the results for the lateral tension (see preceding 
Section), in the same Molecular Dynamics runs. We show S{q) for nominal sizes 
36 X 36, 55 X 55, and 100 x 100. Formula (3.6) represents very well all data obtained 
so far. Formula (3.6) is semiempirical in the sense that the first term results from 
the single-mode approximation to the capillary wave hamiltonian[13,16], whereas 
the polynomial is there to represent the non-singular bulk-like background of <S'(g) 
in the bulk including the liquid solvent. Also we must not forget that at or 
near q = 27r/(T (i.e. near x ~ 40) there appears the nearest-neighbour peak. 
The nearest-neighbour peak is also present in a distorted form in the height- 
height correlation functions. The rise in S{q) towards that peak begins quite 
early, already near a; ~ 4 or even less. The polynomial in (3.6) takes care of that, 
within the limited qf— range of our data. 

The size dependence of S{q) and the lack of it, are illustrated by Fig. 5, Fig. 6, 
and Fig. 7. These are plots of S{q) against x = q^, each for the three sizes at 
equal areas per head a. Fig. 5 is for such a's that 7 ~ 0.84; the systems are clearly 
inside the region (E), as defined in Section II. Fig. 6 is for the tensionless states, 
or almost tensionless, and Fig.7 is for small a, negative 7, thus deep inside region 
(C) of floppy bilayers. 

These three Figures demonstrate the behaviour of the singular term (3.3). As 
reported below, k is always positive and k > 1. Then for positive g as x ^ the 
curve flattens a little as it crosses over from (kx'^ + gx)~^ to l/gx (see Fig. 5). For g 
nearly vanishing (see Fig. 6) this flattening disappears and the rigidity divergence 
1/kx'^ dominates. For A; > 0, ^ < (see Fig.7) there is a pole on the real axis at 

x^x^ = {-g/k) . (3.7) 

Then S{x) diverges faster than it did, aiming at infinity at the asymptote x = x^. 
This is what we see in Fig.7 where each system, depending on its size, has its own 
asymptote. Otherwise at higher values of x, the size appears to affect S{q) very 
little. The poles of (3.3) and (3.6) are discussed further below in subsection B. 
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The data on S{x) were fitted to (3.6). Since one has to be careful with least- 
square fits, as a precaution we have used several versions of (3.6) for each set of 
data, (a) with p ^ and with p = 0, (b) with W2 either zero or not (c) with wi 
zero or not. 

Fits proved to be robust in most cases but gave scattered results for the 
coefficients k and p for systems with large 7 and area a. One firm conclusion is 
that always the coefficient k was positive A; > 0; not a single instance was ever 
found with A; = or A; < 0. Thus k can be interpreted as a rigidity coefficient (up 
to a constant). 

Parameter p was sometimes erratic but other parameters, representing in 
(3.6)the smooth background, did vary smoothly with a and not much. Protrusions, 
accounted for as p/x, are small- wavelength fiuctuations and physically a protrusion 
may be hardly distinguishable from a nearest-neighbour interaction. That may 
explain a correlation of the least-square p with tu^'s. 

Fig 8. shows the variation of the least-square parameter k with the area a. 
Referring to the two regions (C) and (E) described in section II, we can distinguish 
the (E) region of fast increase of k with a, so that the stiff bilayer of high 7 at 
large a has a high rigidity, and the buckling floppy bilayer at low a and negative 
7 in region (C) has a relatively low k. This seems plausible. In region (C) the 
low /c ~ 5. stays constant whereas in region (E) it increases overall, in some fits 
linearly, in some very little. Eliminating a we obtain the variation of least-square 
parameters with 7; Fig. 9 shows the overall increase of k{'y). 

Fig. 10 shows how closely the least-square coefficient g follows 7. Nevertheless 

consistently (7 > 7, suggesting these are two different quantities. In fact, it may 
be argued that g describing the spontaneous distortion-fluctuation of the inter- 
face, is coupled to the hypothetical "true area", whereas 7 is coupled to the "pro- 
jected area" flxed by the experimentalist. This issue was discussed also recently [7]. 
These considerations go back to the distinctions between different molecular ex- 
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pressions[9] for the surface tension of the hquid interface. It appears[17] that in a 
"normal" interface ^ = 7 and in a membrane-hke sheet of a bilayer, g from S{q) 
need not be equal to 7. Consistently, as the Figure shows, ^ > 7. 

B. The presence of a pole on the g— axis. 

Returning now to Fig. 7 which shows S{q) for negative 7 we notice that in- 
creasing the size essentially increases the range as the asymptote x = is 
pushed to the left to lower values. All the divergences are well handled by (3.3) 
and (3.6), but for < 0,A; > the axis is divided into two parts < a; < a;^ 
and x^ < X. For x > x^ all is well, but for < a; < x''" the expressions for S{q) 
become negative and — — 00 as a; tends to from below. Clearly in the inter- 
val < a; < the expression for S is unphysical — as the very definition of the 
scattering factor S makes it a positive quantity. The interval < a; < a;^ must be 
excluded. As a consequence, we cannot consider the limit ^ 0"*" or a; — > 0"*" , as 
is usual for a discussion of the capillary waves and of the "capillary divergence". 

The resolution of these points comes from the realization that in defining S{q) 
as the Fourier transform of a correlation function in the r— space, in fact we are 
dealing with a Fourier series. Indeed the positions r vary continuously but the 
q— vector, {qx,qy), is restricted to 

2tt 

qx = {-r>x (n^ = 0,±l,±2,...) (3.8) 

and similarly for qy. The case = = being excluded, the lowest possible 
value of |g| is go = 27r/Lj; and x cannot be smaller than 

xo = Air^/L^Ly. (3.9) 

Thus the size of the system is ever present in S{q). For positive k and g one can 
speak of the limit g — > 0, but for negative g this is not correct. 

Knowing that for given size x must fulfill a; > a;o we must check that the 
asymptote a; = a;^ is always beyond reach, i.e. that 

0<x^ <xo = Att'^/A (3.10) 
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is always fulfilled. This is the case with all our data and with all our least-square 
fits. Never any least-square fit had to be rejected because the pole would fall in 
the physical region xq < x. 

Now we consider the size dependence of the position of the pole. As the 
parameters k and g vary smoothly (allowing for scatter) so does and therefore 
the calculated values of {—g/k) for some small positive g are also included. 

Of great significance is the plot in Fig. 11 of x^ against 7 because here the 
different sizes fall on a common curve which, with reasonable accuracy, aims at 
the point (0,0). This is very satisfactory: it means that with the increased system 
size the position of the asymptote x = x"^ tends to 0+. Thus the unphysical region 
< x < x^ shrinks to zero with the increased system size. 

Such behaviour appears very satisfactory and also in agreement with the 
content of Section II. 

In this way we have resolved a serious difficulty which was either ignored or 
misinterpreted in the past. 

IV. SUMMARY and DISCUSSION. 

We used atomistic simulations of dimer molecules forming reverse bilayers in 
a solvent. We obtained results on the bilayer isotherm 7(a), the apparent buckling 
transition, the scaling the negative lateral tension if the floppy bilayer, the size 
dependence of the structure factor S{q), and on the new divergence of S. 

These results emphasize again the deep differences between the membrane- 
like two-dimensional sheets of surfactants and the transition region between two 
coexisting fluid phases that is given the name of an interface. 

The bilayer peculiarities were studied in Section II; the lateral tension of a 
stable bilayer though given by (2.1) can be null or negative. This is possible owing 
to the existence of a proper or intrinsic area A of the membrane, constant in a 
flrst approximation. In addition to the thermodynamic parameters, the state of 
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a bilayer depends on the externally enforced area - the edge of the box in the 
simulation. The function 7(a) is the the bilayer isotherm. 

The few simulations that have gone beyond the tensionless state and have 
determined the bilayer isotherm for a large range of arcas[2,4,5,6,7], all found 
interesting and nontrivial behaviour. We find size-independent linear 7 = 26. (o — 
a*) for large 7 > 0. In the Lennard-Jones units (Appendix A) the slope is 26.e/cr^. 
In small systems[l,2,5] the isotherm continues down with little change in slope, 
but in bigger systems there is a very quick change of slope and a transition to 
a floppy state of a buckling bilayer with negative 7 of a flat isotherm saturating 
to the size dependence (2.8). The constraint of constant intrinsic area makes the 
bilayer to buckle. Negative 7's can be extremely close to zero. 

The bilayer together with the solvent fluctuates owing to thermal motion. In 
larger systems the undulations destroy the extension of the linear portion of the 
isotherm into negative 7's; the tensionless state lies in the region of fast change of 
the slope. It is an unwelcome piece of news. It means that the inverse compressib- 
lity Ka taken at the tensionless state, not only is not easily determined accurately 
but moreover is not representative - unless the system is small enough. Rather, 
the slope of the linear part, common to different sizes, is representative. 

We also find that that changeover, from region (E) of the common slope 
to region (C) of negative 7, is very sudden and abrupt, suggestive of a buckling 
transition[18]. However, such a transition, in the strict meaning of a mathematical 
singularity, has not been predicted for a bilayer. We find no definite increase in 
sharpness with size. In bilayers made of long-chain surfactant molecules [6, 7], the 
transition of smaller systems was sudden and appeared truly discontinuous; for 
large bilayers the shape of 7(a) became the same as in Fig.l. A picture emerges of 
a discontinuous transition which is destroyed by the increase in size of the bilayer. 

The structure factor S{q) ~< hqh-q > was independent of size as long as 
a > ao; in the region of floppy bilayer the strong size dependence appears and the 



15 



new pole at q > replaces the capillary wave divergence at q— > 0+. This issue 
is now satisfactorily resolved (see Section IIIB). 

The theoretical picture leading to good predictions [11] for the floppy bilayer, 
is very attractive in its simplicity: as the area A is constant, there is no surface 
tension contribution, only bending. Free energy expressions i.e. mesoscopic hamil- 
tonians are constructed building on such ideas [11]. However it is only but a first 
step; the predicted transition is too smooth, not all size dependence is correctly 
obtained (see Section II), and the structure factor S{q) is not known; it ought to 
appear with another coefficient laxgei than a (see Fig.lO). 

Finally, we remark that the relations between k (i.e. k) and the inverse 
compressibility Ka, taken from the theory of elasticity of plates[l], seem to do 
well in bilayers made of long chains[l,7] but fail for dimers studied here. 
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VI. APPENDIX A . 

The units are those commonly used for systems with Lennard- Jones (6-12) 
potential 

uoo{r)=^e{{a/rY^-{a/rf) {A.l) 

The unit of energy is e and the unit of length is a. Areas ^ or a are given in 
units of (7^, free energy and energy in units of e, reduced temperature T* = kT/e 
where k is the Boltzmann constant, reduced lateral tension 7* = 7cr^/e, reduced 
pressure p* = pa'^/e, compressibility modulus Ka in units of e/a^. In the text all 
quantities are reduced quantities and the asterisk is dropped from the notation. 

In this work all constituent particles are spherical and interact according to 
(A.l). The system contains A^^ "a" solvent particles and A^^ surfactant dimers 
made of A^^ further "a" particles and of Nd "b" particles. The constituents of the 
surfactant "a-b" dimer are bound by an unbreakable but flexible chemical bond. 
Particles are of equal sizes and have equal masses. The potentials are cut and 
shifted[19,20]; for a pair a — P 

Uapir) = iuoo{r) - woo(r^^)) x rjir^^ - r) {A.2) 

where Q!=a,b and /?=a,b. The Heaviside function is ri{x) = 1 or if a; > or 
X <0. For like particles, a = P, the cutoff distance is taken at 

r^^ = 2.5 . (A3) 

(in units of a). For a ^ P, i.e. for "a,b" pairs 

raf3 = ^ 2VV (AA) 

makes the potential purely repulsive in the spirit of WC A [3, 5, 6, 19]. 

The bilayers studied in this work are "reverse" bilayers so named in analogy 
to the reverse micelles. The reverse bilayers were constructed and found to be 
stable[5]; are formed not by entropic effects in conjunction with the hydrophobic 

17 



effect, but primarily by energetic preference. The latter is achieved by making 
the "b-b" attraction much stronger. To achieve this, we modify the energy depth 
parameter e for the pairs "b-b" (only), keeping the common value of e for all other 
pairs. That is 

m = Seaa- 

where the strength parameter S is to be choosen S » 1. Then the "b" -particles, 

the "b"-ends of the surfactant dimcrs, make the strongly cohesive core of the 
bilayer whereas the weakly attracting "a" -ends stay outside in contact with the 
" a" -solvent. 

The Molecular Dynamics simulations were done in a well-established manner: 
at constant volume V, and particle numbers N, N^, by using Verlet leap-frog al- 
gorithm and Nose-Hoover thermostat [19,20]. In all series quoted here the reduced 
temperature was T = 1.9 and the overall density was a reasonably high liquid-like 
density of p = N/V = 0.89204. The augmented "b-b" strength parameter was 
S = 4. The total number of particles contained A''^ dimers that is A''^ "b" parti- 
cles all of them bound in "a-b" dimers, A''^ "a" particles also permanently bound 
in "a-b" dimers, and N — 2Na free "a" particles which were the solvent. Given the 
numbers and the above density p the volume V was fixed; once the area a 

was chosen, A = Naa followed, hence = Ly as the square root, whereas was 
adjusted to fit the volume V = AL^. The areas a, as shown in Fig.l were chosen 
to be in the interval of stable bilayers. 

For the size (iii) (nominally of 36 X 36) 

N = 40000, Nd = 2238, N - 2Nd = 35524, 33. < < 35.5, 41.2 < < 36.6 
For the size (ii) (nominally of area 55 x 55) 

N = 160000, Nd = 5760, N-2Nd = 148480, 51. < < 57.8, 69. < < 53.5 
For the size (i) (nominally of area 100 x 100) 

N = 281344, Nd = 18720, A^ - 2A^d - 243904, 93. < < 105., 36.5 < < 30. 
The proper normalization of the structure factor S{q) leads to the prediction 

18 



< hqh-q >= S{q) = (kT/LxLy) • ^3.2"*^^^ for the singular part of S. The size- 
independent quantity is AS. The Fourier components of h{x,y) were calculated 
for each monolayer separately according to (3.5); the number of molecules in either 
will vary with time. 
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VIII. FIGURE CAPTIONS 



Caption to Fig.l 

The lateral tension 7 plotted against area per dimer head, a, for three sizes : 
size (i) with diamonds, size (ii) with plus signs, size (iii) with squares. All points 
are at T = 1.9, high liquid- like density p = 0.89204; size (i) is that of an area 
about 100. X 100., size (ii) - of an area about 55. x 55., and size (iii) - of an area 
about 36. by 36.. See the Appendix for further details. A is total area, is the 
number of dimer heads, a area per head. Note common curve for large 7 > 0. 

Caption to Fig. 2 

The plateau of 7(a) < at small a, plotted against 1/A = is shown to 
follow the line y ~ 300 * x. See text. The near-constant plateau tends to zero 
with increased system size. Logarithmic plots confirm this Figure. For system 
parameters see Appendix A. 

Caption to Fig. 3 

Total intermolecular energy per particle, U, plotted against a, area per dimer 
head, for the largest system of nominal size 100 x 100. Note a discontinuity, though 
rounded, in the derivative. 

Caption to Fig. 4 

The 7(a) plot for a hypothetical first-order transition in a virtually infinite 
system: fall in 7 with decreasing a towards ao in region (E); a discontinuous change 
of slope to a near- vanishing value in region (C). See text. 

Caption to Fig. 5 

The structure factor S{q) plotted against x = for three system sizes at the 
same area per head a = 1.09 : Diamonds - size (i) 100 x 100; plus signs - size (ii) 
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55 X 55; square boxes - size (iii) 36 x 36. The line shows the fit of the function 
(3.6) to data on system of size (i). The size dependence is seen to be neghgible, 
except for the extension of range of q^. The lateral tensions are strongly positive 
7 = 0.85- 0.90. 

Caption to Fig. 6 

The structure factor S{q) plotted against x = for three system sizes at the 
same area per head a = 1.05 : Diamonds - size (i) 100 x 100; plus signs - size (ii) 
55 X 55; square boxes - size (iii) 36 x 36. The line shows the fit of the function 
(3.6) to data on system of size (i). The size dependence is seen to be negligible, 
except for the extension of range of q^. The lateral tensions are near zero: 7 = 
0.02-0.05. Compare with preceding and following Figure. 

Caption to Fig. 7 

The structure factor S{q) plotted against x = for three system sizes at 
negative 7 and in a floppy state ( area per head a ~ 0.973 — 0.993). Diamonds - 
size (i) 100 x 100; crosses - size (ii) 55 x 55; square boxes - size (iii) 36 x 36. The 
thick line shows the fit of the function (3.6) to data on system of size (i); the two 
other lines are there to guide the eye. The size dependence is seen to be negligible, 
except for the extension of range of q^. The lateral tensions are all negative but in 
the crupling floppy state (C), depend now on size. Compare with preceding two 
Figures. 

Caption to Fig. 8 

The variation of the least-square parameter k with the area a. Note A; > 0. 
Diamonds - size(i) ~ 100, Crosses - (ii), 55, Squares - (iii), 36. Note the two 
regions (E) and (C); (C) buckling floppy bilayer with low and constant k and (E) 
extended tensioned bilayer with k roughly linear with area per head a. 
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Caption to Fig. 9 

The variation of the least-square parameter k with the lateral tension 7. Di- 
amonds - size(i) Lj: ~ 100, Crosses - (ii), 55, Squares - (iii), 36. Note overall 
increase of k with 7. 

Caption to Fig. 10 

The variation of the least-square parameter g with the lateral tension 7. Di- 
amonds - size(i) Lx ~ 100, Crosses - (ii), 55, Squares - (iii), 36. Note how closely 
g follows 7 though systematically g is the larger of the two. 

Caption to Fig. 11 

Pole position plotted against the lateral tension Diamonds - size(i) Lx ~ 
100, Crosses - (ii), 55, Squares - (iii), 36. See text. 

end of Captions, end of paper. 
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